Thermal conductivity of MgO, MgSiC>3 perovskite and post-perovskite in the Earth's 

deep mantle 
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We report lattice thermal conductivities of MgO and MgSi03 in the perovskite and post-perovskite 
structures at conditions of the Earth's lower mantle, obtained from equilibrium molecular dynamics 
simulations. Using an advanced ionic interaction potential, the full conductivity tensor was calcu- 
lated by means of the Green-Kubo method, and the conductivity of MgSiC>3 post-perovskite was 
found to be significantly anisotropic. The thermal conductivities of all three phases were parameter- 
ized as a function of density and temperature. Assuming a Fe-free lower-mantle composition with 
mole fractions a;M g Si0 3 = 0.66 and xm s o = 0.34, the conductivity of the two-phase aggregate was 
calculated along a model geotherm. It was found to vary considerably with depth, rising from 9.5 
W/(mK) at the top of the lower mantle to 20.5 W/(mK) at the top of the thermal boundary layer 
above the core-mantle boundary. Extrapolation of experimental data suggests that at deep-mantle 
conditions, the presence of a realistic amount of iron impurities lowers the thermal conductivity of 
the aggregate by about 50% [2^| . From this result and our thermal conductivity model, we esti- 
mate the heat flux across the core-mantle boundary to be 10.8 TW for a Fe-bearing MgO/MgSi03 
perovskite aggregate and 10.6 TW for a Fe-bearing MgO/MgSiOs post-perovskite aggregate. 



INTRODUCTION 

The thermal conductivity of minerals in the Earth's 
mantle is an important geophysical parameter which gov- 
erns the heat flux from the core up to the surface and 
hence strongly influences mantle dynamics 27 1. More- 
over, the thermal conductivity of minerals at the core- 
mantle boundary (CMB) determines the amount of heat 
extracted from the core, driving the convection of the liq- 
uid outer core and thus controlling the power available to 
the generation of the Earth's magnetic field @, 0|. Yet, 
measuring thermal conductivities at mantle pressures 
and temperatures is extremely challenging, and exper- 
imental data are scarce. Several schemes exist to extrap- 
olate thermal conductivities measured at low er p ressures 
and temperatures to deep-mantle conditions H3] , but 
they are plagued with large uncertainties. Hence a com- 
putational approach is desirable to evaluate thermal con- 
ductivities directly at the relevant conditions. The aim 
of this study is to provide reliable values for the lat- 
tice thermal conductivities of MgO, MgSi03 perovskite 
(Pv) and post-perovskite (PPv) at lower-mantle condi- 
tions and their variation with temperature and density 
(or pressure). These results can be directly applied to 
thermal transport in the lower mantle. 

In deep-mantle minerals, heat is conducted by phonons 
and electromagnetic radiation. The importance of the ra- 
diative contribution to thermal transport in the Earth is 
under debate, and current estimates span a considerable 
range: while [111 ] report a radiative thermal conductivity 
below ~ 0.5 W/(mK) across the lower mantle, [43| pre- 
dict ~ 5 W/(mK) at the CMB, and [H[ even values of up 
to ~ 10 W/(mK), which is of the same order of magni- 
tude as the lattice contribution. Moreover, the radiative 
conductivity seems to depend strongly on crystal grain 



size and on the iron content [13| . In view of these diffi- 
culties, we focus on the lattice contribution in this study. 
If the radiative conductivity turns out to be significant it 
can simply be added to the lattice part presented here. 

Over the past years, different atomic-scale methods 
were developed to calculate lattice thermal conduc- 
tivities. |42| applied the non-equilibrium or "direct" 
method [25j, |28( to derive the thermal conductivity of 
MgO, using molecular dynamics (MD) simulations based 
on density functional theory (DFT). In this approach, 
an energy current from the cold to the hot side of the 
simulation cell is imposed. From this current and the 
steady-state temperature gradient which builds up, the 
thermal conductivity is obtained via Fourier's law. While 
computationally rather efficient, the method suffers from 
strong finite-size effects, thus requiring extrapolation to 
infinite system size and introducing considerable uncer- 
tainties [39(. An approach based on phonon lifetimes, ob- 
tained from DFT, was used by @,0| and by to calcu- 
late the thermal conductivity of MgO. Phonon lifetimes 
were either calculated from line widths in the Fourier 
transform of the velocity autocorrelation function 6] or 
from anharmonic lattice dynamics |47j |. Combined with 
the Boltzmann transport equation for the phonon gas, 
they yield the thermal conductivity in the relaxation time 
approximation. This approach treats the anharmonicity 
of lattice vibrations perturbatively and is thus limited 
to temperatures where atomic displacements from the 
equilibrium positions are small enough for higher-order 
anharmonicity to be neglected. 

A third approach, the Green-Kubo method, uses the 
Green-Kubo relations [2(j to obtain thermal conduc- 
tivities from appropriate current correlation functions, 
which, in turn, are readily extracted from equilibrium 
MD trajectories. This method has been successfully 
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applied to solids (e.g. 0, HH, 48 1) and liquids (e.g. correlation functions in thermodynamic equilibrium: 



It 



33l 38 1 ) . In contrast to the non-equilibrium method, 
no concerns about leaving the linear-response regime 
arise for equilibrium MD. Moreover, the Green-Kubo 
method exhibits a weaker finite-size effect [39|], provides 
the full thermal conductivity tensor in one simulation and 
takes into account thermoelectric effects which can con- 
taminate results of the non-equilibrium method for ionic 
conductors [38j. Unlike the lattice dynamics approach, 
the Green-Kubo method takes into account anharmonic- 
ity to all orders. Thus its validity is not restricted to low 
temperatures. In the light of these advantages, we de- 
cided to use the Green-Kubo approach to calculate ther- 
mal conductivities of MgO, MgSi0 3 Pv and MgSi0 3 PPv 
at conditions spanning a wide pressure and temperature 
range. We also determined conductivities at conditions 
where experimental data are available, and satisfactory 
agreement with these experiments makes us confident 
that our results are equally reliable at CMB conditions. 
A drawback of the method is that it requires long run 
durations (in the nanosecond range) to obtain reasonable 
statistical accuracy. Our calculations are based on clas- 
sical MD simulations involving an interaction potential 
of first-principles accuracy [3] . 
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where ks is Boltzmann's constant, and the are Carte- 
sian components of the energy current (A = U) or of the 
mass currents (A = 1, . . . , N — 1, where N is the num- 
ber of chemical species in the system), with respective 
dimensions of energy or mass times velocity. Angular 
brackets denote an ensemble average. We assume that 
the center of mass is at rest, hence there are only N — 1 
independent mass currents for a system with N chemical 
species. Then, for a system with two species, the thermal 
conductivity is given by (lfl | , 



I 

2^2 



r aa 



aa\2 
Ul) 



T act 



,ae{x,y,z} (3) 



and for a system with three species by [36 
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It is worth noting that equations [3] and 2] are written here 
in terms of mass currents, whereas they were originally 
derived in terms of ionic currents. 



The thermal conductivity tensor A is defined by 
Fourier's law, ]q = — AVT, under the constraint that 
no mass or electric currents are present. This constraint 
is relevant to electronic or ionic conductors, where ther- 
moelectric effects occur [3j. Fourier's law is of linear- 
response type and relates the heat current density jq to 
the temperature gradient VT. For cubic and orthorhom- 
bic crystals, A is diagonal if the coordinate axes are along 
the crystal axes, and direction-dependent conductivities 
can be defined by 



A <* = -Jq/^cT, ae{x,y,z} 



(1) 



In the framework of non-equilibrium thermodynamics 
0, IH , the thermal conductivity can be expressed in terms 
of kinetic coefficients Lab , as is done in equations [3] and 
[4] below. They determine the linear response of the sys- 
tem to deviations from equilibrium, i.e. energy and mass 
flows resulting from thermal and chemical gradients. The 
gist of the Green-Kubo method is that the kinetic coef- 
ficients Lab, although representing non-equilibrium be- 
havior, are linked to fluctuations in thermodynamic equi- 
librium via the fluctuation-dissipation theorem. The ki- 
netic coefficients, and hence the thermal conductivity, 
can therefore be obtained from equilibrium MD by means 
of appropriate Green-Kubo formulae, which relate the 
linear response of a system with volume V to current 



SIMULATION DETAILS 

We performed equilibrium molecular dynamics simu- 
lations in the NVT ensemble, with a time step of I fs 
for the integration of the equation of motion and a Nose- 
Hoover thermostat 



14|, [29| maintaining the system at 



the desired temperature. The cell dimensions were cho- 
sen as the average cell size in a previous NPT run at 
the desired pressure P, maintained by a barostat [24 1. 
The interactions between atoms were described by an 
advanced ionic interaction potential which was param- 
eterized non-empirically, using DFT as a reference [l6j . 
This potential has been shown to reliably predict prop- 
erties of minerals of the system CaO-MgO-Al 2 03-Si02 
over a wide temperature and pressure range, with ac- 
curacy comparable to DFT. In particular, the ionic in- 
teraction potential used in this study was shown to de- 
scribe MgO and the MgSi03 phases perovskite and post- 
perovskite well, predicting lattice constants to within 1% 
and elastic constants mostly to within 10%, compared to 
DFT results [l6| • The elastic constants determine vibra- 
tional modes of the crystal in the limit of long wave- 
lengths 1]. These modes close to the Brillouin zone 
center, in turn, are expected to make the largest con- 
tribution to the thermal conductivity of the crystal 47 1. 
Therefore, we expect the interaction potential to produce 
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accurate lattice dynamics and thermal transport proper- 
ties. For MgO, MgSi0 3 Pv, and MgSi0 3 PPv, we used 
cubic or orthorhombic supcrcells containing 512, 960, and 
720 atoms, respectively. For each composition, tempera- 
ture, and pressure, we generated trajectories of at least 
0.5 ns and up to 2.4 ns. 
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FIG. 1: Thermal conductivity of MgO at 300 K, GPa as a 
function of correlation time r, see eq. [2l averaged over 15 MD 
blocks of 100 ps each. Clearly, a plateau is reached at 30 ps. 

At each time step of the MD run, the mass currents 
for each species and the energy current were extracted 
for later calculation of the different current correlation 
functions needed in eq. [2] An explicit expression for the 
energy current for polarizable ions was derived by [33| . 
The total MD run was then divided into blocks of equal 
length (50 to 100 ps) which were analyzed independently 
for current correlation functions and thermal conductiv- 
ity. The correlation time r in eq. [2] was chosen large 
enough for the thermal conductivity A to reach conver- 
gence. In practice, A as a function of t oscillates around 
its limiting value, and for each MD block, we took a time 
average over the first plateau of the cumulative A(r) (av- 
eraged over the A(r) from the individual MD blocks), see 
Fig. [T] Finally, the thermal conductivity and its ler un- 
certainty were obtained by averaging the results from all 
blocks. 



RESULTS AND DISCUSSION 
MgO 

For reference, we first calculated the thermal conduc- 
tivity of isotopically pure MgO in the fee structure at 300 
K and ambient pressure. At these conditions, the model 
predicts a density p = 3.602 g/cm 3 , which is in excel- 
lent agreement with the experimental density at ambient 



conditions, 3.583 g/cm 3 [4(|. In Fig.[IJ we show the com- 
puted thermal conductivity as a function of correlation 
time r, see eq. The data contain a number of outliers, 
due to the second term on the right side of eq. El Both 
Lm and L\\ are expected to be close to zero when no 
diffusion is present, thus the quotient may take on very 
large positive or negative values occasionally, although it 
should be small in a crystal. Since outliers tend to dis- 
tort arithmetic averages, we calculated the conductivity 
for each MD block from the median of the data over the 
plateau, which is a more representative measure for the 
expectation value in such cases. 

Our value of the thermal conductivity at ambient con- 
ditions is (111 ± 16) W/(mK), significantly larger than 
that found in other computational and experimental 
studies (table J]). However, overestimation with respect 
to experiments is to be expected, since we considered 
a perfect, isotopically pure crystal, whereas the experi- 
ments were performed on real crystals with natural iso- 
topic composition and defects, which reduces the thermal 
conductivity considerably relative to its perfect-crystal 
value [lj| H^. Therefore, our results should indeed be 
larger than the experimental ones. j4?J evaluated the 
isotope effect for MgO and found that at ambient con- 
ditions, the thermal conductivity of an isotopically pure 
crystal exceeds the one of natural samples by as much 
as 46%. This correction for isotopic composition is al- 
ready included in their results in table [I] If we apply the 
same correction to our data, we get A]yigo(300K,0GPa) = 
(76 ± ll)W/(mK). Defects, impurities and grain bound- 
aries in real crystals will further reduce the thermal con- 
ductivity, and thus our result is fully compatible with 
the measured conductivity of (65 ± 15) W/(mK) [l7j |. 
On the other hand, the computed values given by [7| and 
in particular by [42j seem to fall at the low end of values 
reconcilable with experiments, as the computational data 
represent isotopically pure, perfect crystals and therefore 
should not agree with conductivities measured on real 
samples. Following [39j], the relatively small value of [42[ 
may be attributed to the use of a linear extrapolation to 
account for finite-size effects in the non-equilibrium MD 
method, which leads to a systematic underestimation of 
the thermal conductivity. 

The thermal conductivity was evaluated at four more 
p, T points, up to lowermost-mantle conditions, T = 3000 
K, p — 5.307 g/cm 3 , see table |U These data allow us to 
parameterize the behavior of thermal conductivity over a 
wide density and temperature range, including the con- 
ditions relevant to the lower mantle. The temperature 
and density dependence of the thermal conductivity is 
a highly complex matter, and no general theory is cur- 
rently available [32j . The dependence on density can be 
described, in the framework of the Debye approximation, 
as A oc p a , where a is itself a function of density and tem- 

13| tested 



perature in principle 



In a recent study, 



the validity of different models for a(p, T) by measuring 
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TABLE I: Thermal conductivity of fee MgO, MgSiO a Pv, and MgSi0 3 PPv. P M d is the pressure resulting from our MD 
simulations with the indicated density and temperature, and Pe s the ° ne obtained from an equation of state by [4^. Results 
are from this study if not otherwise indicated. For the last p, T point of PPv, the three directional conductivities average to 
(15.3 ± 0.8) W/(mK), to be compared with the calculated bulk A of (15.1 ± 0.9) W/(mK). The difference, which is safely 
within the error bars, results from slightly different averaging schemes: for the bulk value, we first averaged over all directions 
and then over time (i.e. over the plateau, see fig. [1]), whereas for the direction-dependent conductivities, we averaged each 
direction individually over the plateau. 



phase 


II 3\ 

P (g/cm ) 


T (K) Pmd (GPa) P EoS (GPa) 


\ IWT II TSW 

A (W/(mK)) 


MgO 


3.602 


300 





1 


111 ± 16 












76 ± 11 a 












62 












59 ± 6 c 












75 












65 ± 15 c 




5.410 


300 


134 


151 


1400 ± zoO 




4.201 


2000 


41 


45 


40.0 ± 2.5 




5.307 


2000 


133 


148 


141 ± 11 




O.oU 1 


3000 


138 


155 


(D.o =t 4.4 


MgSlC>3 Pv 


4.544 


300 


26 


31 


27.0 ± 2.2 










26 


in f 
19 










31 


10.6 ± 0.6 B 




5.332 


300 


107 


111 


61.3 ± 7.9 










108.4 


23.7 ± 4 g 




4.544 


2000 


40 


42 


9.7 ± 1.0 




5.401 


3000 


139 


137 


12.4 ± 2.0 






300 




ambient 


5.1 h 






300 




ambient 


5.8 i 


MgSiOii PPv 


5.631 


298 


135 


138 


167 ± 25 










141 


65 ± 14 g 




5.482 


2000 


130 


132 


16.8 ± 0.5 




5.631 


2000 


150 


150 


20.6 ± 1.7 




5.482 


3000 


138 


140 


15.1 ± 0.9 



A. r : 
X y : 
A z : 



18.0 ± 
13.7 ± 

14.1 ± 



1.8 
1.0 
1.2 



This study, with isotope correction from [47[ 
4711 . DFT, with isotope correction 
3, DFT, perfect crystal 
7J, DFT, perfect crystal 
17 1, experiment at ambient conditions 
22 1 , experiment at 300 K 

31 1, experiment at 300 K, pressure determined experimentally 

341. experiment at 300 K, metastable (quenched to ambient pressure) 

311 ] . high-P experiments, extrapolated to ambient pressure 



the thermal conductivity of CaGeC>3 perovskite, which 
is an analog phase for MgSiC>3. They found subtle dif- 
ferences in the density dependence of a compared to the 
case of MgO, which they attributed to the larger number 
of optical phonons in the perovskite phase. However, in 
view of the limited number of data points in our present 
study, we take a to be constant, as in [42| . which yields 
an effective a for the entire density range spanned by our 



data points. Concerning the temperature dependence, 
thermal conductivity approximately follows a power law 
A oc T~ b at high temperatures Following Q and [42| . 
we write 

with a fixed reference point po = 3.602 g/cm 3 , Tq = 300 
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density (g/cm ) 



FIG. 2: Density and temperature dependence of the thermal 
conductivity for MgO, Pv and PPv. Data points are results of 
our MD simulations. Lines are fits to eq.[5](MgO, Pv) or eq.|6] 
(PPv), with colors representing different temperatures. Blue: 
300 K, green: 2000 K, red: 3000 K. Inset: PPv at conditions 
relevant to the lower mantle. 



TABLE II: Parameters for the density and temperature de- 
pendence of thermal conductivity resulting from fits to our 
calculated data points, for MgO in fee structure, MgSiOa Pv 
(both described by eq. and PPv, described by eq. [(J] Ref- 
erence density po and temperature To are fixed during the 
fit. 
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MgO 


3.602 


300 


129 


5.42 


1.10 


MgSiOa Pv 


4.544 


300 


28.3 


4.06 


0.608 


MgSiOa PPv 


5.482 


3000 


14.6 


7.48 


0.327 



K, and free parameters Ao, a, and b. The result of a least- 
square fit of eq. [5] to the data points is given in table Ql] 
The quality of the fit can be assessed by means of fig. [5] 
At the high temperatures prevailing at the CMB, an- 
harmonicity, i.e. phonon-phonon scattering, is expected 
to be the dominant mechanism limiting the thermal con- 
ductivity, compared to other sources of phonon scattering 



like isotopic disorder, point defects, and grain bound- 
aries. This is because the number of phonons present 
in the material increases with temperature and hence 
the mean free path between phonon-phonon collisions 
becomes shorter than the mean free path imposed by 
other scattering mechanisms. The reduction of the ther- 
mal conductivity due to isotopic disorder in MgO has 
been shown to decrease from 46% at room temperature 
to only 4% at 4000 K [13]. Qualitatively the same behav- 
ior has been observed in experiments on other materials, 
see e.g. In MD simulations, anharmonicity, i.e. the 
dominant scattering mechanism at high temperatures, is 
automatically included to all orders, and hence our con- 
ductivity results for a perfect crystal should be a good 
approximation for real MgO at CMB conditions. The 
same argument applies to the MgSiOs phases to which 
we turn in the following paragraphs. Although the above 
reasoning is based on sound physical considerations, we 
emphasize that more experimental or simulation data are 
needed to fully understand and quantify the influence of 
different kinds of defects on the thermal conductivity at 
high pressure and temperature. 



MgSiOa perovskite 

MgSiOs in the orthorhombic perovskite structure 
Pbnm is generally accepted to be the most abundant 
mineral in the Earth's lower mantle (i.e. below a depth 
of 670 km). It consists of a three-dimensional network of 
corner-sharing SiOg octahedra, with Mg occupying the 
larger inter-octahedral sites. The calculated thermal con- 
ductivity of MgSiOs in the perovskite structure, averaged 
over all directions, at four state points spanning a wide 
range of densities and temperatures, is given in table U 
along with available experimental data. The effect of iso- 
topic disorder on the thermal conductivity is not known 
for the MgSiOa phases Pv and PPv. As in the case of 
MgO, it may be significant at low temperatures but is 
expected to decrease rapidly with temperature. The den- 
sity and temperature dependence of thermal conductivity 
is well described by eq. [5J and the respective fit parame- 
ters are listed in table|TTl In fig.[H the model conductivity 
is plotted along with the computed data points. 

The available experimental data for perovskite scatter 
considerably (table [I]), and the effect of grain boundaries, 
isotopic disorder and other defects like possible cracks in 
the samples is not known quantitatively. Hence a com- 
parison to our results is difficult. As expected, the results 
obtained from perfect-crystal simulations are larger than 
the experimental values. Our results seem compatible 
with those of 22j but more difficult to reconcile with the 
data of |3l|. Note that at 300 K and approximately 30 
GPa, [31| report a considerably lower conductivity than 
p2j : the data, both derived from measurements, differ by 
a factor of almost 2. At ambient conditions, 34 report 
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the conductivity to be 5.1 W/(mK), which is in line with 
the value of 5.8 W/(mK), derived from a fit to experi- 
mental data by (3l| . 



MgSiC>3 post-perovskite 

The perovskite structure of MgSiC>3 transforms to an 
orthorhombic post-perovskite phase (Cmcm) at approx- 
imately 125 GPa and 2500 K [H, [30] which is believed 
to be stable in the Earth's lowermost mantle close to the 
core-mantle boundary and might be responsible for the 
D" seismic discontinuity (THJl . It is characterized by layers 
of corner- and edge-sharing SiOg octahedra perpendicu- 
lar to the b axis, with Mg occupying inter-layer sites. This 
anisotropic structure exhibits strongly anisotropic elastic 
properties , which should lead to direction-dependent 
phonon velocities, and hence we expect anisotropic ther- 
mal transport properties. Therefore, in addition to the 
direction-averaged conductivity, we also calculated the 
thermal conductivities separately along the three axes of 
the orthorhombic crystal for one data point, at condi- 
tions representative of the lowermost mantle. The simu- 
lation time was extended to 2.35 ns in this case to ensure 
satisfying statistics for each direction individually. The 
results for several p, T values are listed in table HI 

The calculated data points could not be fit satisfacto- 
rily with eq. 03 due to the very weak temperature depen- 
dence of the thermal conductivity above 2000 K (com- 
pare the second and fourth data point of PPv in table [j}. 
This flattening of the thermal conductivity as a function 
of temperature is consistent with experimental observa- 
tions which show a near-constant conductivity above a 
certain temperature, depending on the material studied 
[35| . In parameterizing the thermal conductivity of PPv, 
we therefore used an expression which reconciles a strong 
temperature dependence at lower temperatures with a 
flat behavior at high temperatures. A good fit could be 
obtained with the following functional form: 



A(p,T) = A 



To 



T - 295K 



,; 'Tq - 295K\ b 



To 



V 



(6) 



where po and To are fixed reference values and Ao, a, 
and b are fitting parameters, and the last (normalizing) 
factor ensures that X(po,T ) = X . The results of a least- 
square fit are listed in table Inland compared to the data 
in fig. [21 Due to the somewhat empirical nature of the 
assumed temperature dependence, the validity of eq. |5]is 
restricted to the temperature range covered by our data 
points. It can certainly not be applied below 298 K: in 
fact, the expression diverges at T — 295 K. However, we 
stress that all data points are well fitted. In particular, 
the density and temperature dependence at conditions of 
the lower mantle is well captured by the model, as can 
be seen from the inset in fig. [51 



The calculated conductivity at conditions of the low- 
ermost mantle is clearly anisotropic, and it is lowest in 
the y direction (along the b axis of the crystal). This is 
consistent with the fact that the crystal is softer along b 
(perpendicular to the layers formed by corner- and edge- 
sharing SiC>6 octahedra) than along a and c (in the plane 
of the SiOe sheets). This leads to lower phonon velocities 
along b, at least close to the Brillouin zone center, and a 
reduced conductivity. 

The thermal conductivity of (167 ± 25) W/(mK), ob- 
tained at 300 K and p = 5.631 g/cm 3 , is considerably 
higher than the value derived from experiments at sim- 
ilar conditions by [HI], which is (65 ± 14) W/(mK) (ta- 
bic [l| . Although we cannot quantify the effect of defect 
scattering in their polycrystalline PPv sample of natu- 
ral isotopic composition, the discrepancy seems too large 
to be completely explained by this mechanism, and its 
origin remains unclear. We note however, that in the 
case of perovskite, the thermal conductivity reporte d by 
31 1 was much lower than the experimental value by [22j . 
Interestingly, the estimate of [31| for the thermal con- 
ductivity of PPv at 3000 K and 135 GPa, based on an 
assumed temperature dependence, is close to and even 
slightly higher than our result at similar conditions (19.5 
W/(mK) and 15.1 W/(mK), respectively). This agree- 
ment at high T may partially be due to a fortuitous can- 
cellation of discrepancies, since [3l| assumed a different 
temperature dependence than the one we found. But it 
also hints at the fact that differences between perfect- 
crystal simulations and real-sample experiments become 
less important with increasing temperature. 



IMPLICATIONS FOR THE THERMAL 
CONDUCTIVITY OF THE EARTH'S LOWER 
MANTLE 

When applying our results to heat transport in the 
Earth's deep mantle, they should be considered upper 
estimates for the lattice thermal conductivities of real 
minerals. Our calculations do not take into account the 
natural isotopic composition, impurities, and defects of 



the minerals, all of which lower the conductivity. [22 1 
measured the effect of realistic amounts of iron impuri- 
ties on the thermal conductivity of MgO and MgSi0 3 per- 
ovskite, at relatively low temperatures and pressures. By 
extrapolation, they estimate that 20 mol% and 3 mol% 
of iron in MgO and perovskite, respectively, reduce the 
thermal conductivity of the aggregate at the CMB by 
about 50% relative to that of the chemically pure aggre- 
gate. On the other hand, we did not take into account 
radiative heat transport as an additional mechanism of 
thermal conduction. 

To calculate the thermal conductivity in the lower 
mantle, we assumed a simplified mantle composition, de- 
rived from the pyrolitic composition given by 36], with 
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molar fractions £MgSi0 3 — 0-66 (perovskite structure) 
and XMgO = 0.34. The conductivities were calculated 
along a mantle geotherm, with the depth-dependent pres- 
sure taken from the Preliminary Reference Earth Model 
and the temperature profile adopted from [4l|. Since 
our model for thermal conductivities was parameterized 
as a function of density and temperature, the pressures 
along the geotherm had to be converted to densities. This 
was done by means of equations of state for the mineral 
phases, described by [II], with a revised set of parameters 
for that model taken from (4t|. The pressures resulting 
from the equations of states for a given density and tem- 
perature agree well with the pressures obtained directly 
from our MD simulations (table fl}. This further corrobo- 
rates the adequacy of the interaction potential used in the 
simulations in describing material properties over the /?, T 
range of the lower mantle. By means of the equations of 
state, the mole fractions of the individual mineral phases 
can be converted to volume fractions, yielding about 82% 
Pv and 18% MgO by volume in the lower mantle. These 
numbers change slightly with pressure and temperature, 
and the exact values resulting from the equations of state 
have been used throughout the study. 

The thermal conductivity of a two-phase aggregate, ex- 
pressed in terms of the conductivities, Ai, A2, and volume 
fractions, /1, /2, (/1 + /2 = 1), of the individual phases, 
depends on the geometric details of the assemblage. The 
extreme cases are realized by a structure of alternating 
parallel layers of the two phases, with a heat flux par- 
allel and perpendicular to the layers, respectively. For 
the former case (a "parallel circuit" ) , the conductivity of 
the aggregate is maximum and given by the arithmetic 
or Voigt average A max = f\\\ + /2A2. For the latter (a 
"series circuit"), the conductivity takes on its minimum 
or Reuss average A m i n = (/1/A1 + /2/A2)~ . For other 
geometries, not necessarily built from layers, the conduc- 
tivity of the aggregate will lie within these bounds. 

Fig. [3] shows the thermal conductivity of the 
MgSi03(Pv)-MgO aggregate along the model geotherm 
and, for comparison, of the MgSi03(PPv)-MgO aggre- 
gate close to the CMB. Thermal conductivity increases 
with pressure and decreases with increasing tempera- 
ture, see eq. O With increasing depth, the pressure ef- 
fect dominates over the concomitant temperature rise, 
resulting in a net increase of the thermal conductivity. 
Only close to the CMB, the sharp rise in temperature 
in the thermal boundary layer reverses this trend. At 
2891 km depth, i.e. at the CMB, with P = 136 GPa and 
T = 3739 K, the conductivity of a MgO/MgSi0 3 per- 
ovskite aggregate is predicted to lie between A m j n = 13.7 
W/ (mK) and A max = 19.1 W/ (mK), depending on geom- 
etry, with average A = 16.4 W/(mK). With MgSi0 3 in 
the post-perovskite structure instead, we obtain A m i n = 
13.9 W/(mK), A max = 19.3 W/(mK), and A = 16.6 
W/(mK), i.e. changes are not significant. The value 
for the PPv/MgO aggregate is in good agreement with 
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FIG. 3: Average thermal conductivity A along a model 
geotherm [4lJ for an aggregate with mole fractions ZMgSiOj = 
0.66 (black circles: perovskite structure, red diamonds: post- 
perovskite strucutre) and XM g o = 0.34. The dashed lines 
represent A m in and A max , bracketing the geometry-dependent 
aggregate conductivity (see text). Lines are a guide to the eye. 
The grey-shaded area indicates the range of the Pv/MgO ag- 
gregate thermal conductivity under the assumption that iron 
in the Earth's lower mantle reduces the conductivity of the 
aggregate by 50% 



[3lj ]. who estimate the aggregate conductivity at 4000 
K and 135 GPa to be approximately 16 W/(mK). As- 
suming that a realistic amount of iron impurities reduces 
the aggregate conductivities by 50% [22j, the respective 
average conductivities at the CMB are 8.2 W/(mK) for 
the MgO/perovskite and 8.3 W/(mK) for the MgO/post- 
perovskite aggregate. The influence of iron on the ther- 
mal conductivity is treated approximatively here, and 
more data are needed to better quantify this effect at 
high pressures and temperatures. 

Our results for the thermal conductivity across the 
lower mantle for an iron-free composition are in remark- 
able agreement with (albeit slightly larger than predicted 
by) the thermal conductivity model by [22] which is based 
on an extrapolation of low-P, T experimental data to 
CMB conditions, assuming an aggregate of 20% MgO 
and 80% Pv by volume. On the other hand, the model 
used by to extrapolate available conductivity data to 
high P, T (including iron) yields somewhat smaller val- 
ues, ranging approximately between 4 W/(mK) and 7 
W/(mK) across the lower mantle. Also the estimate of 
the lower-mantle thermal conductivity by [ll[ lies below 
our data: for the lattice conductivity, they assumed an 
iron-free mantle composition (20% MgO and 80% Pv by 
volume) and extrapolated experimental data, finding a 
maximum lattice thermal conductivity A max varying from 
approximately 3 W/ (mK) to 11 W/(mK) across the lower 
mantle. [43| . using the Ross model [37] ] and approxima- 
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tive equations of state to derive high-P, T lattice thermal 
conductivities from available data, also obtained values 
lower than ours, with A max not exceeding 8 W/(mK) for 
an iron-free mantle composition (20% MgO and 80% Pv 
by volume). We emphasize that our data are based on 
simulations directly at lower-mantle conditions and do 
not depend on extrapolations. 

Finally, an estimate for the heat flux across the CMB 
is presented. Given the temperatures T± and T2 at two 
different depths z\ and Z2 as boundary conditions, the 
steady-state heat current density jq is determined by an 
integral form of Fourier's law, 



f 



JQ = 



Z2 - Z\ 



dTX( Pl T) 



'A 



where the two concentric spheres corresponding to Z\ and 
Z2 are locally approximated as parallel planes. Taking 
the temperature at the CMB and 120 km above from 
4l| and neglecting the small density variation across this 
layer, we obtain an average CMB heat flux of 21.5 TW 
for a Pv/MgO aggregate and of 21.2 TW for a PPv/MgO 
aggregate. This estimate is based on a specific thermal 
model of the Earth, and a different temperature profile 
at the CMB would lead to a somewhat different esti- 
mate of the CMB heat flux. Assuming that the pres- 
ence of iron impurities reduces the heat flux by 50% 
[22I ]. it is estimated to be 10.8 TW on average for an 
Fe-bearing Pv/MgO aggregate and 10.6 TW for a Fe- 
bearing PPv/MgO aggregate, with possible variations by 
about ±20%, depending on the geometric details of the 
two-phase assemblage. These values for the CMB heat 
flux are consistent with previous estimates, spanning a 
wide range from 5 TW to 15 TW [H). 



assumes a constant thermal diffusivity \/(pcp) (p: den- 
sitity, cp: specific heat capacity) across the mantle (e.g., 
[46|). which is poorly constrained, moreover. Together 
with the approximation cp — const., this implies A oc p, 
a rather restrictive assumption, to be contrasted with the 
more flexible eq. [5] 

By combining our thermal conductivity results with a 



CONCLUSIONS 



thermal model of the Earth |41[ , the lattice contribution 
to the CMB heat flux is estimated to be about 11 TW 
for a Fe-bearing two-phase aggregate (virtually the same 
with MgSiOa perovskite and post-perovskite). This rel- 
atively high flux is consistent with recent estimates of 
the heat flux required to generate and maintain mantle 
plumes 21] . Due to the large conductivity contrast be- 
tween MgO and the MgSiC>3 phases, the conductivity of 
the two-phase aggregate depends strongly on the aggre- 
gate geometry. Thus, the CMB heat flux may show large 
lateral variations by up to about ±20%. 
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We performed equilibrium MD simulations and used 
the Green-Kubo method to calculate lattice thermal con- 
ductivities of MgO, MgSiOs perovskite, and MgSiOa 
post-perovskite over a wide range of pressure and temper- 
ature conditions relevant to the Earth's deep mantle. To 
our knowledge, these are the first simulation results for 
the MgSiOs phases. Moreover, the thermal conductivity 
of the lowermost mantle has been determined directly, 
without extrapolation from experimental or computa- 
tional low-pressure or low-temperature data and hence 
is free of the inherent uncertainties. 

The data were used to construct a model for thermal 
conductivities as a function of density and temperature, 
which was then applied to the Earth's lower mantle. The 
thermal conductivity was found to increase significantly 
with depth and to decrease steeply across the thermal 
boundary layer above the CMB. These results may be 
used in geodynamic modeling to refine large-scale sim- 
ulations of mantle convection. In this field, one often 
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